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Abstract 

The goal of this paper is to contribute to the field of nondestructive testing by 
eddy currents. We provide a mathematical analysis and a numerical framework for 
simulating the imaging of arbitrarily shaped small volume conductive inclusions from 
electromagnetic induction data. We include in this paper a discussion on data sampling, 
noise reduction, and on probability of detection. 
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1 Introduction 

Nondestructive testing by eddy currents is a technology of choice in the assessment of the 
structural integrity of a variety of materials such as, for instance, aircrafts or metal beams, 
see [To| ]. In this paper we introduce a new eddy current reconstruction method relying on 
the assumption that the objects to be imaged are small. This present study is related to the 
theory of small volume perturbations of Maxwell's equations, see S. It is, however, specific 
to eddy currents and to the particular lengthscalcs relevant to that case. 
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We first note that in the eddy current regime a diffusion equation is used for modeling 
electromagnetic fields. The characteristic length is the skin depth of the conductive object 
to be imaged |lC[ ]. We consider the regime where the skin depth is comparable to the 
characteristic size of the conductive inclusion. 

Using the .E-formulation for the eddy current problem, we first establish energy esti- 
mates. We start from integral representation formulas for the electromagnetic fields arising 
in the presence of a small conductive inclusion to rigorously derive an asymptotic expansion 
for the magnetic part of the field. The effect of the conductive targets on the magnetic 
field measured away from the target is expressed in terms of two polarization tensors: one 
associated with the magnetic contrast (called magnetic polarization tensor) and the second 
with the conductivity of the target (called conductivity polarization tensor). The magnetic 
polarization tensor has been first introduced in [|| in the zero conductivity case while the 
concept of conductivity polarization tensor appears to be new. 

Based on our asymptotic formula we are then able to construct a localization method 
for the conductive inclusion. That method involves a response matrix data. A MUSIC 
(which stands for Multiple Signal Classification) imaging functional is proposed for locating 
the target. It uses the projection of a magnetic dipole located at the search point onto 
the image space of the response matrix. Once the location is found, geometric features of 
the inclusion may be reconstructed using a least-squares method. These geometric features 
together with material parameters (electric conductivity and magnetic permeability) are 
incorporated in the conductivity and magnetic polarization tensors. It is worth emphasizing 
that, as will be shown by our asymptotic expansion, the perturbations in the magnetic field 
due to the presence of the inclusion are complex-valued while the unperturbed field is real. 
As consequence, we only process the imaginary part of the recorded perturbations. Doing 
so, we do not need to know the unperturbed field with an order of accuracy higher than the 
order of magnitude of the perturbation. An approximation of lower order of the unperturbed 
field is enough. 

The so called Hadamard measurement sampling technique is applied in order to reduce 
the impact of noise in measurements. Finally we provide statistical distributions for the 
singular values of the response matrix in the presence of measurement noise and we simulate 
our localization technique on a test example. 

The paper is organized as follows. Section 2 is devoted to a variational formulation of 
the eddy current equations. Section 3 contains the main contributions of this paper. It 
provides a rigorous derivation of the effect of a small conductive inclusion on the magnetic 
field measured away from the inclusion. Section 4 applies a MUSIC-type localization to 
the eddy current model. Section 5 discusses the effect of noise on the inclusion detection 
and proposes a detection test based on the significant eigenvalues of the response matrix. 
Section 6 illustrates numerically on test examples our main findings in this paper. A few 
concluding remarks are given in the last section. 

2 Eddy Current Equations 

Suppose that there is an electromagnetic inclusion in R 3 of the form B a = z + aB, where 
B C R 3 is a bounded, smooth domain containing the origin. Let T and T a denote the 
boundary of B and B a . Let denote the magnetic permeability of the free space. Let 
fi* and cr* denote the permeability and the conductivity of the inclusion which are also 
assumed to be constant. Wc introduce the piecewise constant magnetic permeability and 
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electric conductivity 



yu* in B a , 
fio in B c a - 



cr a {x) = 



cr» in B a , 
in B%. 



Let (E a , H a ) denote the eddy current fields in the presence of the electromagnetic inclu- 
sion B a and a source current Jo located outside the inclusion. Moreover, we suppose that 
Jo has a compact support and is divergence free: V • Jo = in R 3 . The fields E a and H a 
are the solutions of the following eddy current equations: 



VxE a = iLj/j, a H a in R 3 , 
VxiJ Q = a a E a + J in R 3 , 
E a (x)=0(\x\- 1 ) 1 H a {x) = 0{\x\- x ) 



(2.1) 



as x 



Eliminate H a in we obtain the following ^-formulation of the eddy current problem 



Vx^'Vx^a — \Luo a E a = iujJq in R 3 , 

\7-E a = in B c a , 

E a (x) = 0(|a;| _1 ) as \x\ -> oo. 



(2.2) 



Throughout this paper, let denote the limit values of u(x ± tn) as t \ 0, where n is the 
outward normal to T a , if they exist. We will use the function spaces 



u 



= < u 



e i 2 (R 3 ), Vxm € L 2 (R 3 ), V-m = in B c a 



and 



m : it e X Q (R 3 ) 



n = 



and the sesquilinear form on X Q (R 3 ) x X Q (R 3 ) 

a a {E,v) = (/i^VxE.VxuJgs - iw<7*(.E, u) Ba , 

where (-, ■)£> stands for the L 2 inner product on the domain D C R 3 . The weak formulation 
of the .E-formulation (EDI) is: Find E a g X Q (R 3 ) such that 



a a {E a ,v) =kj(J ,v) B o a , V«eX a (R 3 ). (2.3) 
The uniqueness and existence of solution of the problem (|2.3|) is known (cf., e.g., Ammari 



ct al. Q and Hiptmair 16 ). Note that the constraint L u + ■ n = in X Q (R 3 ) only serves 
to enforce the uniqueness of E a in B c a jl6]]. This is not essential for the validity of the 
^-formulation of the eddy current model. We have 



a a {E a ,v) = iu(J ,v) B c, V«eX a (R 3 ). 



(2.4) 
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Throughout the paper we denote by Eq the unique solution of the problem 

Vx/i 'VxEo = iojJ inR 3 , 

V-£ = in R 3 , (2.5) 
E (x) = 0(|a;| _1 ) as \x\ -> oo. 
The field E satisfies 

^VxEo.Vxd^ =iu(J ,v) v , Vw £ H^cur^IR 3 ), (2.6) 

where H_i(curl;R 3 ) = (it : = £ L 2 (R 3 ) 3 ,Vxm 6 L 2 (K 3 ) 3 ). 

3 Derivation of the Asymptotic Formulas 

In this section we will derive the asymptotic formula for H a when the inclusion is small. 
Let k = w/xocr*. We are interested in the asymptotic regime when a — > and 

v = ka 2 (3.1) 

is of order one. Moreover, we assume that fx* and no are of the same order. 

In eddy current testing the wave equation is converted into the diffusion equation, where 
the characteristic length is the skin depth <5, given by i5 = yj2jk. Hence, in the regime 
v = O(l), the skin depth S is of order of the characteristic size a of the inclusion. 

We will always denote by C a generic constant which depends possibly on fx^/no, the 
upper bound of uj^Loa^a 2 , the domain B, but is independent of w, cr*, ^o, /•**■ Let \i r = /Lt*//io- 

3.1 Energy Estimates 

We start with the following estimate. 

Lemma 3.1 There exists a constant C such that 

\\Vx(E a - Eo)\\ L »p») + ^WEa - Eu\\ L 2 {Ba) < Ca 3/2 (Vk\\E \\ L ™(B a ) + ||Vx£^,|| l » (Bq) ). 



Proof. By (gj) and (|2J), we know that 

(fj,~ 1 Vx(E a - E Q ),Vxv) M 3 - iu}(a a (E a - E ),v) Ba 

= Guq 1 - n- 1 )(WxE ,Wxv) Ba +iuj(<T a E ,v) Ba , VveX a {R 3 ). (3.2) 



Since 
and 



\(VxE ,Vxv) Ba \ < Ca 3 / 2 ||Vx^||iao (Bo) ||Vxi;|| £a( B a 
\{a a E ,v)\ < Ca 3 / 2 a4E \\ L ^ {Ba) \\v\\ L 2 {Ba)7 



by taking v = E a — Eo £ X Q (R 3 ) in (3.2) and multiplying the obtained equation by //o we 
have that 

/i^HVx^ - E )\\ 2 L2m + k\\E a - Eo\\l HBa) 
< Co?' 2 (\^xE a \\ L ^ {Ba) \\Vx(E a - E )\\ L 2 {Ba) + k\\E \\ L oo (Ba) \\E a - E Q \\ L 2 (Ba) ) . 
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This completes the proof. 



Let H 1 (B a ) = {ip G L 2 {B a ),Vip G L 2 {B a f}. Let <j) G H x (B a ) be the solution of the 
problem 

- A^o = —V ■ F in B aj -d n <f>a = (Efo(as) - F(x)) • n on T Q , / <j) dx = 0, (3.3) 
where 

F[x) = i[V 2 x x (x - z) + i[D z (V z x E )(z)}(x - z) x (x - z). (3.4) 

Here [D Z (V Z x Eo)(z)]ij = d Zi (V z X Eo(z))j is the (i, j)-th element of the gradient matrix 
D z (S7 z xE )(z) ofV z xE (z). Let tr denote the trace. Since tr[D(Vx#o)] =V'(Vx£ ) = 
0, we know that 

V x F(x) = V z x E (z) + [D Z (V Z x E Q )(z)](x - z). (3.5) 
Note that since E$ is smooth in B a we have 

||V x E Q - V x F\\ L ~ (Ba) < Ca 2 \\V x E \\ w ^ {Ba) . (3.6) 

Denote by Hq = (kj/io)~ 1 Vx.Eo an d introduce w G X Q (R 3 ) as the solution of the 
problem 

a a (w,v) = iw^o(/^o 1 _ ^^(Hoiz) + DH {z)(x - z), Vx«) Bq 

+ioj(a a F,v) Ba , Vd £ X Q . (3.7) 
The following lemma provides a higher-order correction of the error estimate in Lemma [O. 



Lemma 3.2 Let w be defined by (3^7). There exists a constant C such that 

\\Vx(E a -Eq- w)\\ L 2 {R s } < Ca 7 ' 2 {\l - ^\ + z/)||V x E \\ W 2^ {Ba} , (3.8) 
\\E a -E Q - V^o - to|U»(B ) < Ca 9/2 (|1 ~ AV 'I + f)||V x ^|| wa ,. (B<>)l (3.9) 



where v is given by ( p.l[ ). 

Proof. First we set ■i/' G iJ 1 (B Q ) and 17 be such that g = oiiT a , Ag = in and 
9 = O(|o;| _1 ) at infinity. Let 

J W> inB Q , 
" '~ 1 V 5 in 



Since v G X Q , it follows from (2.4) that 



iw(a aJ E a ,VV)B =0, VV> G H\B a ). 
This yields V-E a = in B a and 2£~ • n = on T a . 
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Similarly, we can show from (3.7) that w ■ n = —F{x) ■ n on r Q and V-w = — V • F in 
B a . From (|]|) we also know that V-(S +Vc/> ) = V-FmB a and (S o +V0 o )~ n = F{x)-n 
on r a . Thus 

V • (S a - E - V0 O -w)=0 in B a , (E a - S - V0 O - w)~ ■ n = on r Q , 

which implies by scaling argument and the embedding theorem that 

\\E a - E - V0 O - w\\ L 2 {Ba} < Ca\\V x (E a - E -Vcpo-w)\\ L 2 {Ba) 

= Ca\\Vx{E a -E -w)\\ L 2 (Ba) , 



for some constant C independent of a and <t*. Therefore, (3.9) follows from ( p.8| 
To show (|3.8|), we define 0o as the solution of the exterior problem 



— A(po = in .B^, 0o = </>o on r Q , <po — > as |x| — > oo. 

The existence of O in W l ~ l {B c a ) = \cp : = G L 2 {B c a ),VLp G L 2 (B^) 3 } is known 

(cf, e.g., Nedelec ||). 

Define # = V0 O in B«, *o = V0 O in then $ G X a (R 3 ). It follows from dO) and 



(13.71) that 

(^~ 1 Vx( J E Q - S - #0 - w),Vxo)i3 - iw(cr Q (S Q - S - * - w),v) Ba 
= iw/i (Mo 1 - /i^X-ffo - fo(«) - Dflo(2)(i - 2),Vxi)) Bt> + iu;(er Q (S + *o - f),w) Ba . 
Multiply the above equation by fiQ we have then 

(H0^a 1 ^ /x ( E a -E0-&0- W), VX») B S - ifc(S Q - S - *0 ~ W,v) Ba 

= - /V X )(Ho - Ho(*) - DH (z)(i - z),Vx))) Bo + ifc(S + * - F, v) Ba .(3.10) 

It is easy to check that 

\iwfi (H Q - H Q (z) - BH (z)(x - z),Vxv) Ba \ < Ca 7/2 \\iujnoH Q \\ W 2, x . {Ba) \\V xv\\ L 2 {Ba) 

= Ca 7/2 \\V x Eo\\ W 2,oo {Ba) \\Vxv\\ L2(Ba) . 

Now taking v = E a - S - * - w G X Q (R 3 ) in ( |3.10| ), since V x $ = in I 3 and 
<J?o = V</>o in B ai we obtain that 



|Vx(S a -S -™)" 2 



k\\E a -E -V4> ~wf 



| L 2 (R 3) T n,||£/ Q - ^0 - v TO - <"IIl 2 (B ) 

< Ca 7 / 2 |1 - n~ 1 |||V X JSo|| W ».«(B a ,)||VXi;||La(B < ,) + fc||S - S + V^ollz^i^lMU 2 ^) 

< Ca 7 / 2 (|1 - /i- 1 ! + !/)||V X ^llwa-CB^IIVx^lU^fl,,). 
Here, we have used 

1 1 So - F + X7<p \\LHB a ) < Ca\\S7 x (S - F)|| £a(Ba) < Ca 9 / 2 ||V x S ||„^ (Bo) (3.11) 

and ||i>||l2( Sq ) < Ca||V x v\\ L 2( B \, since E — F + V^o and v are divergence free in B a 
and have vanishing normal traces on r a . This shows (3.8) and completes the proof. □ 
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We note that T)Hq(z) is symmetric since V x Hq(z) = 0. Hence, by Green's formula, 
Oo 1 - nZ x )(H {z) + T>H Q (z){x - z),Vxv) Ba 

= Oo 1 - /C 1 ) J t (( H o(z) + T>H (z){x - z)) x n) ■ vdx 

= J [fi- 1 {H {z)+-DH {z){x-z))xn]r a -vdx, 

where [-]r a stands for the jump of the function across T. Let w(£) = w(z + a£), we know 
from that, Vu £ Xi(R 3 ), 

(/x _1 Vxt/), Vx«)|3 — iuja 2 (aw, v)b = iacJ/Uo J {H (z) + aT>Ho(z)£) x n]r ■ vd£ 

+iuja 2 (aF(z + a£),v) B , 

where fx(£) = ^ if £ £ B, fi{£) = /i if ^ G B c and a(£) = cr* if $ e B, a{£) = if £ £ B c . 
This motivates us to introduce the solution wq(£) of the interface problem 

x /i _1 Vg x Wo — \Loaa 2 w = itoaa 2 \a~ x F{z + a£)] in B U B c , 

V 4 • wo = in S c , 

[too x n]r = 0, [/i _1 V^ x w x n]r = — iw(l - fi^ 1 )(H Q (z) + aDH (z)£) x n on L, 
I «>o(€) = Od^l" 1 ) as |£| oo. 

/ 3? — J£ 

It is easy to check that w(x) = awo I 

\ a 

The following theorem which is the main result of this section now follows directly from 



(3.12) 



Lemma 3.2 



Theorem 3.1 There exists a constant C such that 

< Ca 7/2 (\1 - Hr 'I + ")I|V x E \\ W 2,oe (Ba) , 



x — z 

Vx [E a - E - aw ( ) 

a 



L 2 (l 



E a - E - V0 O - aw 



a 



L' 2 (B a 



<Cc^' 2 {\l- ^~ 1 \ + v)\\V x E, 



Oil W 2 <°°(B a )- 



To conclude this section wc remark that 



a _1 F(z + a£) = iu^o [ - H {z) x £ + -T>H (z)£ x £ 



= ic^o \\ Yl Ho ^ ei x f £ Dffo(*)«eiej€ x £ ] (3.13) 

where DH'o(- 2 )ij is the (i, j/ )-th element of the matrix DUqCz) and T denotes the transpose. 
Thus 



w (£) = iufto Uj2 H ^^ + f £ Dfro(«)< 3 -*«(€) 

\ *=1 i,j=l y 



(3.14) 



r 



where is the solution of the following interface problem 

x /i _1 V$ x Oi \Loaa 2 9i = iujocPei x | in B U B c , 
V € • 0i = in 5 C , 

[0i x n] r = 0, [^ _1 V^ x OiX n] r = — [^i _1 ] r e i x n on T, 
I e i (0 = O(|^|- 1 )as|e|^oo, 

and V&ij (£) is the solution of 

x /i _1 V^ x — \Luoa 2 ^ij = iuj<7a 2 ^jei X^inBU B c , 
V € • = in B c , 

[i&ij x n] r = 0, x ty^ x n]r = — 3[/t _1 ]r£jej x n on I\ 

[ *. tJ (0 = 0(|^|- 1 )as|^|^ O o. 



(3.15) 



(3.16) 



Here is the unit vector in the Xi direction. It is worth emphasizing that since v = O(l), 
6i and are uniformly bounded in Xi(R 3 ). 

We impose V • Oi = outside B to make the solution 9i unique outside B. In this case 
by H, Proposition 3.1] we can show that Oi = 0(|£|~ 2 ) and V X 0i = 0(|£|~ 3 ) as l£l ~> °°> 
which implies by integrating ( 3.15 ) over B that 



; / (6i +e,x |)<f£ = nx ^ _1 V x 

ifl J dB 

= / n x /i _1 V x Oid$, 

JdB R 



— > as i? — > +oo, 
where Br is a ball of radius R so that i? C Br. Thus we obtain 



/ {O i + e i x£)d£ = Q. 
Jb 



(3.17) 



Similarly, by imposing V-\&.y = outside B we know that V x \E^ 7 - = Q (|$|~ 3 )- Moreover, 
integrating ( |3.16[ ) over i? and using similar argument leading to (3.17) together with the 
symmetry of T)Hq(z) yields 



3 

V DH (*)« / + fce< x ejdg = 0. 



(3.18) 



3.2 Integral Representation Formulas 

The integral representation is similar to the Stratton-Chu formula for time-harmonic Maxwell 
equations (cf., e.g., Nedelec H). 



<s 



Lemma 3.3 Let D be a bounded domain in R 3 with Lipschitz boundary whose unit outer 
normal is n. For any E G H„i(curl;R 3 \L>) satisfying VxVxB = 0,V--E = in R 3 \D, 
we have, for any x G M. 3 \D, 

E(x) = -V.X/ (E(y) x n)G{x 1 y)dy - [ (V y x E{y) x n)G{x,y)dy 

-V„ / (E(y)-n)G(x,y)dy, 

where G(x, y) = 4v \x-y\ is the fundamental solution of the Laplace equation. 



Proof. For the sake of completeness we give a sketch of proof. Since E G H_ i (curl; 
for any F such that F(y) = 0(|y| _1 ) and T)F(y) = 0(|y|~ 2 ) as \y\ — > oo, we can obtain 
by integrating by parts, the conditions Vx Vx£ = 0,V'£ = 0in R 3 \5, that 

(E, —AF) M 3\d = (£,VxVxF-VV-f) K3 \ S 

= - (E x n) -VxFdy - \7xExn-Fdy + (E ■ n)\7-Fdy. 

JTo JTd JTjj 

Now for x G R 3 \L> and j G {1,2,3}, we choose F(y) = G(x,y)ej and thus —A y F = 
S(x,y)ej, where S(x, ■) is the Dirac mass at x. Then we have 

Ej(x) = - (E(y) x n) • V y x (G{x 1 y)e j )dy -I (V y x E(y) x n) J G(x,y)dy 
+ [ (E { y).n) d -^l d y 

= - { V x x / (E{y) x n)G(x,y)dy) -[ (V y x E(y) x n) j G{x,y)dy 

[ (E(y)-n)G(x,y)dy, 

where we have used the fact that 

(E(y) x n) ■ V x x (G{x,y) ej ) = -(V« x {G(x,y)E(y) x n)) r 
This completes the proof. □ 

The following lemma will be useful in deriving the asymptotic formula in next subsection. 
Recall that H a = VxE a ,H = — Vx£ . 

Lemma 3.4 Let H a = H a — Ho. Then we have, for x G B%, 

H a {x)= V x G(x,y)xV y xH a (y)dy + (l- f —) (H a (y) -V ' y )V x G(x,y) dy. 
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Proof. It is easy to check that VxiT a = and V-H a = in B c a . By the representation 



formula in Lemma 3.3 we have 



H a (x) = -V. x jf x n)G(x, y)dy -V,| • n)G(x, y)dy, 

where H£ = H a \Bg- Denote H~ = H a \B a and let be defined likewise. By the interface 
condition [E a x n]r a = 0, we have 

H+ ■ n = — !— V x E+ ■ n - H ■ n = — !— div r (-E+ x n) - H ■ n 

M* TT I T 

= — H a n- H n. 

Mo 

where divr a denotes the surface divergence. Then since [H a x n]r a = 0, we have 

H a (x) = -V x x J (H~(y) x ri)G(x,y)dy 

-V. / {^H-{y)-n-H Q {y)-n)G{x,y)dy. (3.19) 
Jr Q Mo 

For the first term, 

-V x x y (H~(y) x n)G{x,y)dy 

x (H a (y)G(x,y))dy 

J B a 

= V x / (G(as, y)V y x H^y) + V y G(a:, y) x H a {y))dy 

= j B (y x G(x,y)xV v xH a (y) + (H a -V !C )V y G(x,y))dy, (3.20) 

where we have used the identity 

Vx (u x ») = ti(V-u) — (it • V)u + (w • V)it — u(V-u), 

and the fact that V x • W y G(x,y) = —A y G(x 1 y) = 0. For the second term, we first notice 
that 

-V. / (—H~(y) ■ n - H Q (y) ■ n)G(x, y)dy = -^vj H~(y) • nG(x,y)dy 
Jv a Mo Mo Jr Q 

+ (1 _ £i) Va! f H (y) ■ nG(x, y)dy. 
Mo Jv a 



10 



By integration by parts we have 

V x / H a (y) ■ nG(x,y)dy = W x I V y • (G(x,y)H a (y))dy 



Similarly 



Thus 



V* / V y G(x, y) ■ H a (y) + G(x, y)V-H a (y)dy 



(H a (y)-W y )W x G(x,y)dy. 



V, / (H (y) ■ n)G(x, y)dy = [ (H Q {y) ■ S7 y )S7 x G{x 7 y)dy. 

-V x f (^H-(y)-n-H (y)-n)G(x,y)dy 

Jr a Mo 



Mo Jb 



(H a (y)-V y )V x G{x,y)dy 



+ (1 - — ) / (Ho(y) ■ V v )V x G{x, y)dy. 

MO JB a 

This completes the proof by substituting ( 3.20| )-(3.21) into (3.19). 



(3.21) 



3.3 Asymptotic Formulas 

In this subsection we prove the following theorem which is the main result of this section. 

Theorem 3.2 Let v be of order one and let a be small. For x away from the location z of 
the inclusion, we have 



H a (x)-H (x) 



3 

\Y,H a (z)i f DlG( X ,z)^x(di + eiX^ 

A 1 *r B 



Wfi-^ 

M* 



H (z) q J D 2 x G(x, z)(ei + x 6^ 



R(x), 



where (D x G)ij = dl iXj G and 

\R(x)\<Ca 4 \\H Q \\ W 2^ iBa 
uniformly in x in any compact set away from z. 
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Proof. The proof starts from the integral representation formula in Lemma 3.4. We 



first consider the first term in the integral representation in Lemma 3.4. By Theorem 3.1 
we know that 



\E a -E - V0 O - av)o(?—Z)\\ L » (Ba) < Ca 9 / 2 (\l - /z" 1 ] + i/)||V x E \\ w ^ (Ba) . (3.22) 



, x — z , 
a 

Since VxHq = and VxH a = aE a in B a , we have 



V x G(x,y) x V y x H a (y)dy = a* / V x G(x,y) x E a (y)dy = h 



where 



y- * 



Ii = cr»y V x G(a;,y) x ^ Q (y) - #o(y) - V0 o (y) -atw ( 
h = v* [ V x G(:e, y) x (^(y) + Vfo(l/) - F(v))dy, 

JB a 

Is = a* J b (V*G(x,y)-V x G(x,z)-DlG(x,z)(y-z))x (F(y) + aw (l^)yy, 



U = <r* / (V x G(a;,z)+£>iG(aj,2:)(y-z)) x F(y) + ati» ( 



y - z 



dy. 



By ( |3.22D , we have 



III < Ca 6 (|l-A*r 1 |+"V«l|Vx^|| W 2,oo (B(>) 

< G^ll-^'lll-ffolk^^) 

< Ga 4 j|i?o|| w 2 <=(B„)- 



By (|T3h we have |I 2 | < Ga 6 er*||V x £ || W 2,«, (Sc( ) < Ca 4 ||i? ||vi/^oo (Bci) . Similarly, by 
using (3.4) and (3.14) we can show |I 3 | < Ca 4 ||i?olliv 2 '° c (B„)- F° r the remaining term we 
first observe that 

I 4 = ia 4 (7* f (V x G(x,z) + aDlG(x,z)£) x (a^Fiz + a£) + w Q (£))d£. 
On the other hand, 

1 3 



which implies after using (3.18) 

3 

h 



ika 4 



+ika 5 



-Y^Hoiz), / W x G(x,z) x (e. t x £ + 9i)d£ 
1 i=i Js 

i ^ H (*)< ^ r£G(x, z)i x ( e , x € + 0,)^ 



fli(x), 
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where |J?i(a;)| < Ca 4L \\H \\w 2 -< x (B a )- Using (3.17), this shows that 



V x G(x,y) x V y x H a (y)dy 



ika 5 



3 

\ £ H„(z) s I D 2 x G(x, z)t- x (e, x £ + fc)^ 



+ H 2 (x), (3.23) 



where 



|i* 2 (x)| < Cfca 6 |l - ^HlHolk^CB,,) < Ca A \\H \\ w ^ {Ba 



Now we turn to the second term in Lemma |3.4j. From Theorem 3.1| we know that 

a 



H a — —H 



Let 



Then 



where 



x — z 
-V x x w ( ) 



L 2 (B a 
1 



< C« 7 / 2 (|l - A*" 1 ) + ^)||^o|| 
V c x w (£). 



/ (H a -V y )V x G(x,y)dy = - f DlG(x,y)H a (y)dy = Hi 



II 4 



Hi = - / DlG(x,y)(H a (y)-^H (y) 

JB a V M* 

H 2 = / (DSG( a!> y)-D2( !B) *))(Ho(y) + HJ(^-5))d V) 

M* Jb " 



Ha 



Mo 
M* 
Mo 



/ DlG(x,z)(H (y)-H (z))dy, 
Jb„ 



,y-z. 



II 4 = / ££G(*,z)(ff (*) + Ho*( i: -^ I ))dV- 

M* Jb q " 



It is easy to see from ( 3.24| ) that |IIi| < Ca 4: \\Ho\\w 1 -°'(B a )- By (3.14) we know that 
\\H*{y^)\\ L , (Ba) <Ca^\\H a \\, 



\W 2 *°°(B a ), 



which implies |II 2 | < Ca 4 \\H \\ W 2 .^(B a )- Similarly, we have | II3 1 < Ca 4 ||Ho||w 1 .°=(.B a )- 



Finally, by (3.14), we have 

3 . 

II 4 = -^VVffo(z)* / DlG(x,z)(e i + -Vxe i )d£ + R 3 (x), 
M* ~[ Jb \ 2 J 

where I-R3 (a?) | < Ca 4 ||ff || w 2 >°°(B a )- Therefore, 
/ (H a -W y )W x G(x,y)dy 

JB„ 



-^VV) -»<,(*)< f DlG(x,z)(e i + lvxe i )d£ + R4(x) (3.25) 
M* l=1 Jb \ l / 
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with | J?4 (as) | < Ca 4 \\Ho\\ W 2.< X ( Ba y This completes the proof by substituting (3.25) and 
(3.23) into the integral representation formula in Lemma |3.4|. □ 



Let M be the 3x3 matrix whose column vectors are | J B £ x (9i+ei x£)d£, | J B £ x (9 2 + 
62 x £)d£, and \ J B £ x (63 + 63 x We call M the conductivity polarization tensor. The 

tensor whose column vectors are (l-^ 2 -) J B ( e i + ^V xO^Jd£, / B (e 2 +|Vx0 2 j^, 

and ^1 — J B ^3 + x O^j d£ is called the magnetic polarization tensor. It reduces in 
the zero conductivity case (a — 0) to the one first introduced in [0. 

Assume that (j,q = fJ-*- The asymptotic formula derived in Theorem 3.2 reduces in this 
case to the following result. 



Corollary 3.1 Assume that fj,Q = fJ>*. Then we have 

H a (x)-H (x) = ika 5 DlG(x,z)M T H (z) + R{x). 



(3.26) 



The remainder satisfies |J2(aj)| < Ca ||i?o||w 2 ' o °(B Q ) uniformly in x in any compact set 
away from z . 

Furthermore, we can show the following symmetry property. 
Proposition 3.1 Assume that (Iq — /x*. Then the tensor M is symmetric. 
Proof. From 

€ x (8i + tn x £) • e, = (8i + tn x £) • {e, x £), 

we obtain that 

/ £ x (6i + e< x £) • e,d£ = — / V x V x t ■ (e 3 ■ x £)d£ 



and therefore, 

£ x (0j + e, x £) • e 3 d£ 
Integrating by parts yields 
I x (0i + e l x £) ■ ejd£ 



1 



V x V x Oj ■ V x V x 0i - ivOj • V x V x Gid£ 



1 

i;/ 



VxVx0j'VxVx0i- ii/V x 0j ■ V x 0;d£ 
+ y (V x Bi x re) - • (6j)-di. 
Since [Vx8,xn] r = and [0j x re] = 0, we have 

y (V x 0, x n)- ■ {6j)-d£ = J (V x 0, x re)+ • (0j) + d£, 
and since V x V x 0.; = in B c , we find that 

/ (V x 0, x re)" • {6j)-d£ = - / Vxfl.-Vx 0^. 

Jr JB" 
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It follows that 



[ £ x (0t + e 4 x |) • e 3 -d^ = ^- / vxVx^-VxVx^- ii/V x 3 - ■ V x 0jd£ 
is w Ub 

V x 6>. t ■ V x 3 -d£, 



which shows the symmetry of M. □ 

Now we assume that Jo is a dipole source whose position is denoted by s 

J Q (x) = Vx(pS(x,s)), (3.27) 
where S(-, s) is the Dirac mass at s and the unit vector p is the direction of the magnetic 



dipole. The existence and uniqueness of a solution to (2.1) follows from j23J. In the absence 



of any inclusion, the magnetic field Ho due to Jo(x) is given by 

H (x) = VxVx(pG(x,s)) = D 2 x G{x,s)p. (3.28) 
We note that Jo is not in the dual of X Q (R 3 ) , however we can form the difference E a — Eq 



and solve for that difference in X Q (IR 3 ). That way we are able to recover Theorems 3.1 and 
The asymptotic formula ( 3.26| ) can be rewritten as 



q (H a - H )(x) ~ ika 5 (D 2 x G(x, z)q) T M (D 2 x G(z, s)p) . (3.29) 



where M is defined in Corollary |3.1[ 



It is worth emphasizing that Proposition 3.1 shows that the reciprocity property for 



(3.30) holds. If the dipole Jo is located at x, then the field p ■ H a at s is the same as the 
one obtained if Jo is located at s and p ■ H a is measured at x. 
Next, writing 

M = SfteM + i3m M, 

we obtain 

5Re(q ■ (H a - H )(xj) ~ -ka 5 (DlG(x, z)q) T (3m M) (DlG{z, s)p) , 

and 

3m(q • (H a - H )(x)) ~ ka 5 (D 2 x G(x, z)q) T (§?e M) (D 2 x G(z,s)p) . 
Since Ho is real, it follows that 

$Sm(q-H a (x)) ~ ka 5 (D 2 x G{x, z)q) T (KeM) (D 2 x G(z, s)p) . (3.30) 



Formula (3.30) will be used in the section for locating and detecting the target. 



4 Localization and Characterization 

In this section we consider that there are M sources and N receivers. The mth source is 
located at s rn and it generates the magnetic dipole Jg"^(r) = V x (p5(r,s rn )). The nth 
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receiver is located at r n and it records the magnetic field in the q direction. The (n, m)-th 

entry of the N x M response matrix A = (A»»)n=i,...,ff,ra=i m is the signal recorded by 

the nth receiver when the mth source is emitting: 



The response matrix is the sum of the unperturbed field Hq (r n ) ■ q and the perturbation 



is the sum of the unj 

ffi™'(r n ) ■ q — H^ n \r n ) ■ q. This perturbation contains information about the inclusion 
but it is much smaller (of order a 3 ) than the unperturbed field. Consequently, it seems that 
we need to know the unperturbed field with great accuracy in order to be able to extract 
the perturbation and to process it. In practice, such an accuracy may not be accessible. 
However, we know that the unperturbed field is real while the perturbation is complex- 
valued, as shown by ( 3.30| ). The imaginary part of the response matrix is therefore equal to 



the imaginary part of the perturbation and this is the data that we will process: 

A = CAo,nm)n=l,...,JV,m=l,...,M, A).nm = 3m (H^ m) (r n ) ■ q) . (4.1) 

We assume that N > M, that is, there are more receivers than sources. As in ||, in 
order to locate the conductive inclusion z + aB we can use the MUSIC imaging functional. 



We focus on formula ( |3.30| ) and define the MUSIC imaging functional for a search point z s 

by 



2mu(-z"" 



1 



Eli II 0-n - P) (DgG(n,zS)q ■ e h . . . , D*G(r N , X S)q ■ e,) T || 2 



1/2 

, (4-2) 



where P is the orthogonal projection on the range of the matrix A and (ej., 62,63) is an 
orthonormal basis of M 3 . From it follows that the following proposition holds. 

Proposition 4.1 In the presence of an inclusion located at z, the matrix Aq has three sig- 
nificant singular values. Moreover, the MUSIC imaging functional Zmu attains its maximum 
approximately at z s = z. 

Once the inclusion is located we can compute by a least-square method the tensor 5fte M 
associated with the inclusion from the response matrix Ao. Given the location of the in- 
clusion, we minimize the discrepancy between the computed and the measured response 
matrices. For a single frequency, knowing -Re M may not be sufficient to separate the con- 
ductivity of an inclusion from its size. However, SfeM obtained for different frequencies w 
may be used to reconstruct both the conductivity and the size of the target. This will be 
the subject of a forthcoming work. 



5 Noisy Measurements 

In this section we consider that there are M sources and N receivers. The measures are 
noisy, which means that the magnetic field measured by a receiver is corrupted by an additive 
noise that can be described in terms of a real Gaussian random variable with mean zero and 
variance o~\. The recorded noises are independent from each other. 
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5.1 Hadamard Technique 

Standard acquisition. In the standard acquisition scheme, the response matrix is mea- 
sured at each step of M consecutive experiments. In the mth experience, m = 1, . . . , M, the 
mth source (located at s m ) generates the magnetic dipole Jq (r) = V x (pS(r,s m )) and 
the N receivers (located at r„, n = 1, . . . , N) record the magnetic field in the q direction 
which means that they measure 

Aneas.nm = ^0,nm + W nm , U = 1,...,N, m = 1, . . . , M, 

which gives the matrix 

A mcas = A + W, (5.1) 



where Ao is the unperturbed response matrix (4J) and W nm are independent Gaussian 
random variables with mean zero and variance <j\. Here, ifa™' 1 (r„) is the magnetic field 
generated by a magnetic dipole at s m and measured at the receiver r„ in the presence of 
the inclusion. 

The so called Hadamard noise reduction technique is valid in the presence of additive 
noise and uses the structure of Hadamard matrices. 

Definition 5.1 A Hadamard matrix H of order M is a M x M matrix whose elements are 
— 1 or +1 and such that H T H = MIm- Here Ijy/ is the M x M identity matrix. 

Hadamard matrices do not exist for all M . A necessary condition for the existence is that 
M = 1,2 or a multiple of 4. A sufficient condition is that M is a power of two. Explicit 
examples are known for all M multiple of 4 up to M = 664 [^2| . 

Hadamard acquisition. In the Hadamard acquisition scheme, the response matrix is 
measured during a sequence of M experiments. In the mth experience, m = 1, . . . , M, all 
sources generate magnetic dipoles, the m' source generating H mm ij'^ n \r). This means 
that we use all sources to their maximal emission capacity (which is a physical constraint) 
with a specific coding of their signs. The N receivers record the magnetic field in the q 
direction, which means that they measure 



M 

H mm ,A 0inm/ +W nm = (A H T ) nm +W nm , n = l,...,N, m = l,...,M, 

m' — l 

which gives the matrix 

B mcas = A H T + W, 

where Ao is the unperturbed response matrix and W nm arc independent Gaussian random 
variables with mean zero and variance o~\. The measured response matrix A moas is obtained 
by right multiplying the matrix B moas by the matrix jjH: 

A mcas = — B moas H = — A H T H + — WH, 

which gives 

A mcas = Ao + W, W = ^WH. (5.2) 



17 



The benefit of using Hadamard's technique lies in the fact that the new noise matrix W has 
independent entries with Gaussian statistics, mean zero, and variance a 2 /M: 

~ ~ 1 M a 2 M 

^[WnmWn'm'] = 2_j ^qmH q ' m 'E[W nq W n r q/ ] = H qm H q , m ,8 nn ,8 qq , 

q,q' = l g,5f'=l 

cr 2 M a 2 

^ I j ^ It',!;,'- II )m'q8mi' — \ j - ' ^ ^ H) m ' m (5 nn ' 
9=1 

£_iA A 

where E stands for the expectation and S mn is the Kronecker symbol. This gain of a factor 
M in the signal-to-noise ratio is called the Hadamard advantage. 

5.2 Singular Values of a Noisy Matrix 

We consider in this subsection the case where there is measurement noise but no inclusion 
is present. The measured response matrix is the N x M matrix 

A mcas = W, (5.3) 

where W consists of independent noise coefficients with mean zero and variance (J 2 /M and 
the number of receivers is larger than the number of sources N > M. This is the case when 
the response matrix is acquired with the Hadamard technique and there is no inclusion in 
the medium. 

We denote by a^* > > cr^ 1 ^ > ■ ■ ■ > a^' the singular values of the response 

matrix A moas sorted by decreasing order and by A^ M ) the corresponding integrated density 
of states defined by 

A (M) ([a, b}) = ^Card j/ = 1, . . . , M , a[ M) <= [a, 6]| , for any a < b. 

The density A^ M ) is a counting measure which consists of a sum of Dirac masses: 

1 M 

M ^ 3 

3=1 

For large N and M denote 

7 := N/M > 1. 

The following statements hold. 

Proposition 5.1 a) The random measure A^ M ) converges almost surely to the determin- 
istic absolutely continuous measure A with compact support: 

A([a u ,a v ])= / — p 1 [ — jda, < a u < a v , (5.4) 



where p 7 is the deformed quarter- circle law given by 

— J ((7 1 / 2 + l) 2 - <J 2 ) {(J 2 - (7 X /2 - 1)2) ^ 7 l/2_ 1<(r < l/2 + lj 
7TO~ * 



p 7 (<r) 



otherwise. 



(5.5) 
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b) The normalized I 2 -norm of the singular values satisfies 
1 M 

— y^(c| A/ ' ) ) 2 — 70"n ^y^r/a^Z in distribution, (5-6) 



A/ 

Ml 

where Z follows a Gaussian distribution with mean zero and variance one. 



c) The maximal singular value satisfies 



- r 7 1/2 + 1 + ^273 (! + 7~ 1/2 ) 1/3 ^i + o( jj^Ta )] in distribution, (5.7) 
where Z\ follows a type-1 Tracy-Widom distribution. 
The type-1 Tracy-Widom distribution has the pdf prwi : 

P(Zi < z) = I pTwi(x)dx = exp ( — — / </?0e) + — z)( y 9 2 (x)c?a; S ) , 

where </? is the solution of the Painleve equation 

ip"(x) = xip(x) + 2ip(x) 3 , <p(x) x ^~°° Ai(x), (5.8) 

Ai being the Airy function. The expectation of Z\ is E[Zi] ~ —1.21 and its variance is 
Var(Zi) ~ 1.61. 

Proof. Point a) is Marcenko-Pastur result |18|| . Point b) follows from the expression of 
the normalized / 2 -norm of the singular values in terms of the entries of the matrix: 



M N M 

J_ Vfo- (JVf)> ) 2 - — trf A T A 1 - — W A 2 

j^j / A 7 / — Jtf \ meas-^measy — / J / _ n meas,nm' 

j — 1 n—1 rn—1 



and from the application of the central limit theorem in the regime M 3> 1 ■ The third point 
follows from Jl7||. □ 



5.3 Singular Values of the Unperturbed Response Matrix 



We now turn to the case where there is one conductive inclusion in the medium and no 
measurement noise. The measured response matrix is then the N X M matrix Ao defined 

by 

A , nm = ka 5 (D 2 x G(r n , z)q) T VteM (D 2 x G(z, s m )q). (5.9) 

Let us introduce the singular value decomposition of the symmetric 3x3 matrix 5fte M: 

5ieM=(V M ) T £ M V M , 

where S M is a diagonal matrix with singular values > ajf > af 1 > 0. The matrix Ao 
possesses three nonzero singular values given by 

M 



Ao 



te 5 of[^|(V M (^G(^ Sm ) 9 )) 

m— 1 

TV 

x[Y,\(V™(DlG(r n ,z)q)) 



n=l 



1/2 



1/2 



J = 1,2,3. 
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5.4 Singular Values of the Perturbed Response Matrix 

The measured response matrix using the Hadamard technique in the presence of an inclusion 
and in the presence of measurement noise is 

A mcas = A + W, (5.10) 

where Ao is given by ( |5.9D and W has independent random entries with Gaussian statistics, 
mean zero and variance a 2 /M. 

We consider the critical regime in which the singular values of the unperturbed matrix 
are of the same order as the singular values of the noise, that is to say, a l °, the first singular 
value of Ao, is of the same order of magnitude as cr n . We will say a few words about the 
cases much larger or much smaller than cr n after the analysis of the critical regime. 



Proposition 5.2 a) The normalized I 2 -norm of the singular values satisfies 

M 

M [— ^{<jf f) ) 2 - 70*] AI ^° (a£°) 2 + y/2yo*Z in distribution, (5.11) 
i=i 

where Z follows a Gaussian distribution with mean zero and variance one and 

CT o A ° = [E«°) 2 ] 1/2 - (5-12) 
j'=i 

bl) If erf' < 7 1 / 4 er n , then the maximal singular value satisfies 

a[ M) -^[Z 2 + l + _l 7 ^(l + 1 -y^ Zl +o{^j l )\ m distribution, (5.13) 

where Z\ follows a type-1 Tracy Widom distribution. 
b2) Ifof a > 7 1/4 <J n , then 

a[ M) = of (a + 0{—^)j in probability, (5.14) 

where 

/ 2 4 \ 1/2 

a= (l + (l+7)-xSr+7-^T) • (5-15) 

//, additionally, > , then the maximal singular value in the regime M ^> 1 

has Gaussian distribution with the mean and variance given by 

E[aH = ^(a + oij^)), (5.16) 

Var(a< M) ) = g(/? + (l)), (5.17) 



where 



1-7 



0\4 



1/2' 



(5.18) 
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Proof. Point a) follows again from the explicit expression of the Z 2 -norm of the singular 
values in terms of the entries of the matrix. Point b) in the case N = M is addressed in ]l3| 
and the extension to N > M is only technical and can be obtained from [|l2| [l5|. □ 

Note that, in the item b2), if erf~° = > , then the fluctuations are not Gaussian 



anymore, but they can be characterized as shown in |L3|. Note also that formula (5.18) seems 
to predict that the variance of the maximal singular value cancels when o~^~° \ 7 1 / cr n , but 
this is true only to the order M _1 , and in fact it becomes of order Af~ 4 / 3 . Following 

J we can anticipate that there 
4 er n + !»M -1 ' 3 for some fixed w. 



11 we can anticipate that there are interpolating distributions which appear when o~ 1 ° 



5.5 Detection Test 

The objective in this subsection is to design a detection method which comes with an 
estimate of the level of confidence, in the presence of noise, in our ability to determine 
whether there actually is a conductive inclusion. 

Since we know that the presence of an inclusion is characterized by the existence of three 
significant singular values for Ao , we propose to use a test of the form R > r for the alarm 
corresponding to the presence of a conductive inclusion. Here R is the quantity obtained 
from the measured response matrix defined by 



M-3(l+ 7 - 1 /2)2 Z^j=4 



(M) 
1 

EAi I 1 



R =~ 1 ^172' ( 5 - 19 ) 



and the threshold value r has to be chosen by the user. This choice follows from Neyman- 
Pearson theory as we explain below. It requires the knowledge of the statistical distribution 
of R which we give in the following proposition. 

Proposition 5.3 In the asymptotic regime M»l the following statements hold. 

a) In absence of a conductive inclusion (Eq. $5J{ )) we have 

* - 1 + 7- 1/2 + ^m^ 1/2 (1 + 7- 1/2 ) 1/3 % + o(^), (5.20) 
where Z\ follows a type-1 Tracy Widom distribution. 

b) In presence of a conductive inclusion (Eq. ( 5.1(\ )) : 
bl) If erf' > 7 1 / 4 er n ,, then we have 

erf /3 1 / 2 

where Zq follows a Gaussian distribution with mean zero and variance one. 
b2) I/crf- < 7 1/4 cr n , then we have \5.2<\ ). 
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Proof. In absence of a conductive inclusion, we have on the one hand that the truncated 
normalized P-norm of the singular values satisfies 



M 



M 

5Z( <T j M ' ) ) 2 — 7 a *> ^y^y/j^z m distribution, 



M -3(1 + 7 -V 2 ) 2 



where Z follows a Gaussian distribution with mean zero and variance one, which implies 
that 



1 



A I 



M - 3(1 + 7- 1 / 2 ) 2 



1/2 , 1 



) 2 = 7 V2 ^n + o(^) in probability, (5.22) 



and on the other hand the maximal singular value satisfies (5/7) in distribution. Using 
Slutsky's theorem, we find the first item of the proposition. 

In presence of a conductive inclusion, we have on t he one hand that the truncated nor- 
malized ? 2 -norm of the singular values sat isfie s ( |5.22 ). On the other hand the maximal 



singular value is described by Proposition ^2j which gives the desired result by Slutsky's 
theorem. □ 

The data (i.e. the measured response matrix) gives the value of the ratio R. We propose 
to use a test of the form R > r for the alarm corresponding to the presence of a conductive 
inclusion. The quality of this test can be quantified by two coefficients: 

- The false alarm rate (FAR) is the probability to sound the alarm while there is no inclusion: 

FAR = P(R > r\ no inclusion). 

- The probability of detection (POD) is the probability to sound the alarm when there is an 
inclusion: 

POD = ¥{R > r\ inclusion). 

It is not possible to find a test that minimizes the FAR and maximizes the POD. However, 
by the Neyman-Pearson lemma, the decision rule of sounding the alarm if and only if R > r$ 
maximizes the POD for a given FAR 5 with the threshold 

rs = 1 + 7 - 1/2 + ^rj 27 3 7- 1/2 (l + 7 - 1/a ) 1/a *^(l - (5-23) 

where $twi is the cumulative distribution function of the Tracy- Widom distribution of 
type 1. The computation of the threshold r$ is easy since it depends only on the number of 
sensors N and M and on the FAR 6. Note that we should use a Tracy- Widom distribution 

table. We have, for instance, $twi( - 9 ) - °- 45 ' $ twi( - 95 ) - °- 98 and $ twi( - 99 ) - 2 - 02 - 
The POD of this optimal test (optimal amongst all tests with the FAR 5) depends on 
the value a l and on the noise level a n . Here we find that the POD is 



POD = $ VM 



1/2 

-a — 7 ' r$ 



£1/2 

where $ is the cumulative distribution function of the normal distribution with mean zero 
and variance one. The theoretical test performance improves very rapidly with M once 
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erf" 1 > 7 1/4 cr„. This result is indeed valid as long as erf > 7 1/4 cr n . When erf < 7 1/4 cr n , 
so that the inclusion is buried in noise (more exactly, the singular values corresponding to 
the inclusion are buried into the deformed quarter-circle distribution of the other singular 
values), then we have POD = 1 — $twi(^twi(1 — ^)) = ^- Therefore the probability of 
detection is given by 

a -i 

i i 

-a — 7 ' r$ 



POD = max|$^VM ^ j,SV. (5.24) 

The transition region erf" ~ 7 1 / 4 er n is only qualitatively characterized by our analysis, 
as it would require a detailed study of the statistics of the maximal singular value when 
erf = 7 1/4 cr n + aiA/" 1 ' 3 for some fixed w. 

Finally, the following remark is in order. The previous results were obtained by an 
asymptotic analysis assuming that M is large and a 1 and cr n are of the same order. In the 
case in which <j 1 is much larger than er n , then the proposed test has a POD of 100%. In 
the case in which a l is much smaller than cr n , then it is not possible to detect the inclusion 
from the singular values of the response matrix and the proposed test has a POD equal to 
the FAR (as shown above, this is the case as soon as erf < 7 1 / 4 cr n ). 



6 Numerical Experiments 

In this section, we will give some numerical examples to illustrate the performance of the 
detection algorithm. The unperturbed measurement is acquired synthetically by asymptotic 



formula (3.26) and noisy measurements are given by (5.1C). Assume that B a is an ellipsoid 



described by 



/ \2 , / \2 i ( z z o) 2 2 
{x - x ) +{y- yo) + = a , 
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where a is characteristic length of the inclusion measured in meters. Then the domain B is 
characterized by letting a = 1 and (xo, yo, zq) to be origin. We assume that the inclusion B a 
is also located at the origin, a = 0.01, /Lt* = fj,o = 1.2566 x 10~ 6 H/m and a = 5.96 x 10 7 S/m. 
We let oj = 133.5 to make ka 2 = 1. We compute the conductivity polarization tensor M by 
an edge element code for ( |3.15| ). The numerically computed tensor M is given by 

'2.6185 + 0.350H 0.0000 + O.OOOOi 0.0001 + O.OOOOi \ 
0.0000 + O.OOOOi 2.6180 + 0.3500i -0.0001 - O.OOOOi J . 
^0.0001 + O.OOOOi -0.0001 - O.OOOOi 1.6403 + 0.1930i / 

We remark here that the polarization tensor M is computed numerically and is a diagonal 
matrix for an ellipsoidal inclusion. 

The configuration of the detection system includes coincident transmitter and receiver 
arrays uniformly distributed on the square [—2,2] x [—2,2] x {1}, both consisting of 256 
(M = N = 16 2 ) vertical dipoles (p — q = e^) emitting or receiving with unit amplitude. 
The search domain is a box [— 0.5,0. 5] 3 below the arrays. It is worth mentioning here that 
the number of transducers should be a multiple of 4 in order to be able to implement the 
Hadamard technique in a realistic situation. 

In the above setting, we calculate the SVD of the unperturbed response matrix Ao. 
Figure |l| displays the logarithmic scale plot of the singular values of Aq. Wc observe that 
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our numerical results agree with our previous theoretical analysis: there are three significant 
singular values associated with the inclusion. Then we can construct the projection P with 
the first three singular vectors corresponding to the first three significant singular values. 
In the right part of Figure we also plot the magnitude of Imu on the cross section z = 0, 
which shows that the MUSIC algorithm can detect the inclusion with high resolution. 

We test the influence of the noisy measurements by adding a Gaussian noisy matrix 
with mean zero and variance cr^/M to unperturbed response matrix Ao- In our tests, the 
Gaussian noise is generated by MATLAB function randn. The imaging results shown in 
Figure ^ indicate that the i magin g results become sharper as the noise level is smaller. 
Then we show the validity of ( 5.24 ). Noticing that M = N makes 7 = 1 in our setting. By 
the analysis in Section 5, for given FAR 5, POD depends on the ratio a x /a n . Here we only 
consider the critical regime in which erf -0 is of the same order of u n (specially erf" > u n ). 
Fixing FAR <5, for each ratio a 1 °/cr n , we generate 1000 Gaussian noisy matrices with mean 
zero and variance o\)M and add them to Ao to get according noisy response matrices A. 
We compute R with the help of SVD for each A and count the times for R > r$ to get 
the numerical POD. Figure || shows the comparisons between numerical POD and ( 5.24 ) for 
each 8. We can conclude that the numerical results arc in good agreement with (5.24). 



7 Concluding Remarks 

In this paper we have provided an asymptotic expansion for the perturbations of the mag- 
netic field due to the presence of a small conductive inclusion. This was done under the 
assumption that the characteristic size of the inclusion is of the same order of magnitude as 
the skin depth. Our asymptotic expansion is valid for arbitrary target shapes. That asymp- 
totic formula was in turn used to construct a method for localizing conductive targets. We 
also presented numerical simulations for illustration. Thinking ahead, it appears that it 
would be very interesting to apply the findings from this paper to real-time target identi- 
fication in eddy current imaging using the so called dictionary matching method. We are 
also interested in investigating target tracking from induction data at multiple frequencies. 
In the presence of noise, another problem of interest is to study how to estimate resolution 
for the localization of targets. This will be the subject of a forthcoming publication. 
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